##########################################
########################################## 
## STUDY 5
##########################################
########################################## 
setwd(getwd())
source("libraries.R")

# load data for study
data_raw <- haven::read_dta(here("Data", "study5.dta"))

##########################################
########################################## 
## Recode Variables
##########################################
########################################## 
df <- data_raw %>% 
  mutate(female = case_when(
    Q1 == 1 ~ 0,
    Q1 == 2 ~ 1,
    TRUE ~ NA_real_
  )) %>% 
  mutate(white = case_when(
    Q2 == 2 ~ 1,
    TRUE ~ 0
  )) %>% 
  rename(edu = Q3) %>% 
  
  #Dummy indicating whether participants were in "control" or "social exclusion" condition
  mutate(exclusion = case_when(
    Q7_intro == 1 ~ 0,
    Q8_intro == 1 ~ 1,
    TRUE ~ NA_real_
  )) %>% 
  
  #Need for Chaos
  mutate_at(vars(Q9_1:Q9_8), ~na_if(., 6)) %>% 
  mutate(nfc = rowMeans(select(., Q9_1:Q9_7), 
                        na.rm = TRUE)) %>% 
  
  #Status-Driven Risk-Taking
  mutate_at(vars(Q5_63, Q5_69, Q6_1, Q6_5, Q6_7), ~flip(.)) %>% 
  mutate(sdrt = rowMeans(select(., Q5_63:Q6_7), 
                         na.rm = TRUE))  %>% 
  #Political Rumors: Believing and Sharing of DEMocratic and REPublican rumors
  mutate(true_dem = rowMeans(select(., c("Q10_1", "Q11_1", "Q12_1")), na.rm = TRUE)) %>% 
  mutate(share_dem = rowMeans(select(., c("Q10_2", "Q11_2", "Q12_2")), na.rm = TRUE)) %>% 
  mutate(true_rep = rowMeans(select(., c("Q13_1", "Q14_1", "Q15_1")), na.rm = TRUE)) %>% 
  mutate(share_rep = rowMeans(select(., c("Q13_2", "Q14_2", "Q15_2")), na.rm = TRUE)) %>% 
  
  #Standardize variables
  mutate(across(.cols = c(edu, nfc,sdrt, true_dem, share_dem, true_rep, share_rep),
                ~ (scale(.) %>%  as.vector),
                .names = "{.col}_z_score")) %>% 
  #Scale variables to range from 0 to 1
  mutate(across(.cols = c(edu, nfc,sdrt, true_dem, share_dem, true_rep, share_rep),
                ~ (zero1(.) %>% as.vector),
                .names = "{.col}01"))

##########################################
########################################## 
## Sample Description 
## (See Supplementary Materials SM1, Table S2)
##########################################
########################################## 
#Gender
print(paste(round(mean(df$female , na.rm = TRUE), 2)*100 , "percent of participants are female"))
#Education
print(paste("Less than high school:", round(prop.table(table(df$edu)), 2)[1]*100, "percent"))
print(paste("High school graduate:", round(prop.table(table(df$edu)), 2)[2]*100, "percent"))
print(paste(" Some college, but no degree:", round(prop.table(table(df$edu)), 2)[3]*100, "percent"))
print(paste("2 year college degree:", round(prop.table(table(df$edu)), 2)[4]*100, "percent"))
print(paste("4 year college degree:", round(prop.table(table(df$edu)), 2)[5]*100, "percent"))
print(paste("Graduate or professional degree:", round(prop.table(table(df$edu)), 2)[6]*100, "percent"))
#Ethnicity
print(paste("White/caucasian:", round(prop.table(table(df$ Q2)), 2)[2]*100, "percent"))

##########################################
########################################## 
## Factor Analysis, etc.: Need for Chaos 
## (See Supplementary Materials SM2, Table S7)
##########################################
########################################## 
cfa <- df %>% 
  select(Q9_1:Q9_7) %>% 
  na.omit()

m1  <- ' nfc_oneDim  =~ Q9_1 + Q9_2 + Q9_3 + Q9_4 + Q9_5 + Q9_6 + Q9_7  
        Q9_3~~Q9_4
        Q9_3~~Q9_5
        Q9_4~~Q9_5'

m1 <- lavaan::cfa(m1, data=cfa) 
m1_fit <- summary(m1 , fit.measures=TRUE) 

print(paste("Need for Chaos - Old,  Chi-square:" , round(m1_fit$FIT["chisq"],2 )))
print(paste("Need for Chaos - Old, CFI:" , round(m1_fit$FIT["cfi"],2 )))
print(paste("Need for Chaos - Old, TLI:" , round(m1_fit$FIT["tli"],2 )))
print(paste("Need for Chaos - Old, RMSEA:" , round(m1_fit$FIT["rmsea"],2 )))

##########################################
########################################## 
## Remove 5% with highest Need for Chaos? (SM8)
## When set to TRUE, used for producing Fig S8, Fig S9
##########################################
##########################################
# Set to TRUE for excluding top 5% of participants with highest Need for Chaos
exclude_obs <- FALSE

if (exclude_obs == TRUE) {
  
  df <- df %>% 
    filter(nfc_z_score < quantile(nfc_z_score, probs = .95, na.rm = TRUE))
  
  print("Excluding top 5% observations with highest levels of the Need for Chaos")
  
}

##########################################
########################################## 
## Analysis - Fig. 3 of Main Text
##########################################
########################################## 
df_reduced <- df %>% 
  select(share_dem01, share_rep01, true_dem01, true_rep01, nfc_z_score, female, edu01, white, exclusion) %>% 
  filter(exclusion==0) %>% 
  na.omit()

m_share_dem <- lm(share_dem01 ~ nfc_z_score + female + edu01 + white, data = df_reduced )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "nfc_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study5")

m_share_rep <- lm(share_rep01 ~ nfc_z_score + female + edu01 + white , data = df_reduced )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "nfc_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study5")

m_true_dem <- lm(true_dem01 ~ nfc_z_score + female + edu01 + white , data = df_reduced )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "nfc_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study5")

m_true_rep <- lm(true_rep01 ~ nfc_z_score   +  female + edu01 + white , data = df_reduced )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "nfc_z_score" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study5")

study5_coef <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study5_coef,file=here("Data_recoded","study5_fig3.R"))

#regression table for supplemental materials (Table S12)
stargazer(m_share_dem , m_share_rep , m_true_dem , m_true_rep , 
          align=TRUE,
          no.space=TRUE,
          column.labels=c("Share - Dem", "Share - Rep" , "True - Dem" , "True - Rep"),
          covariate.labels = c("Need for Chaos (std.)" , "Women" , "Education" , "Age  (std.)" , "Ethnicity (1 = white)" ,  "Intercept"),
          keep.stat = c("n","ser"),
          dep.var.labels.include = FALSE)

##########################################
########################################## 
## Analysis - SM5: Social Marginalization Experiment (Table S21)
##########################################
########################################## 
reg1 <- lm(nfc01 ~ exclusion, data = df)
reg2 <- lm(nfc01 ~ exclusion + sdrt01 + female + white + edu01, data = df)

stargazer::stargazer(reg1 , reg2,
                     no.space=TRUE,
                     column.labels=c("NFC" , "NFC" ),
                     keep=c("exclusion" , "sdrt01" , "female" ,"white", "edu01" , "Constant" ),
                     covariate.labels = c("Social Exclusion (1 = Yes)", "Status-Driven Risk-Taking" , "Gender (1 = Women)" , "Ethnicity (1 = white)", "Education" , 
                                          "Intercept"),
                     keep.stat = c("n","ser"),
                     dep.var.labels.include = FALSE
)






